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Abstract 

Background: Recent studies on the medical treatment of Parkinson's disease (PD) led to the introduction of the so 
called Deep Brain Stimulation (DBS) technique. This particular therapy allows to contrast actively the pathological 
activity of various Deep Brain structures, responsible for the well known PD symptoms. This technique, frequently 
joined to dopaminergic drugs administration, replaces the surgical interventions implemented to contrast the 
activity of specific brain nuclei, called Basal Ganglia (BG). This clinical protocol gave the possibility to analyse and 
inspect signals measured from the electrodes implanted into the deep brain regions. The analysis of these signals 
led to the possibility to study the PD as a specific case of dynamical synchronization in biological neural networks, 
with the advantage to apply the theoretical analysis developed in such scientific field to find efficient treatments to 
face with this important disease. Experimental results in fact show that the PD neurological diseases are 
characterized by a pathological signal synchronization in BG. Parkinsonian tremor, for example, is ascribed to be 
caused by neuron populations of the Thalamic and Striatal structures that undergo an abnormal synchronization. 
On the contrary, in normal conditions, the activity of the same neuron populations do not appear to be correlated 
and synchronized. 

Results: To study in details the effect of the stimulation signal on a pathological neural medium, efficient models 
of these neural structures were built, which are able to show, without any external input, the intrinsic properties of 
a pathological neural tissue, mimicking the BG synchronized dynamics. 

We start considering a model already introduced in the literature to investigate the effects of electrical stimulation 
on pathologically synchronized clusters of neurons. This model used Morris Lecar type neurons. This neuron model, 
although having a high level of biological plausibility, requires a large computational effort to simulate large scale 
networks. For this reason we considered a reduced order model, the Izhikevich one, which is computationally 
much lighter. The comparison between neural lattices built using both neuron models provided comparable 
results, both without traditional stimulation and in presence of all the stimulation protocols. This was a first result 
toward the study and simulation of the large scale neural networks involved in pathological dynamics. 
Using the reduced order model an inspection on the activity of two neural lattices was also carried out at the aim 
to analyze how the stimulation in one area could affect the dynamics in another area, like the usual medical 
treatment protocols require. 

The study of population dynamics that was carried out allowed us to investigate, through simulations, the positive 
effects of the stimulation signals in terms of desynchronization of the neural dynamics. 

Conclusions: The results obtained constitute a significant added value to the analysis of synchronization and 
desynchronization effects due to neural stimulation. This work gives the opportunity to more efficiently study the 
effect of stimulation in large scale yet computationally efficient neural networks. Results were compared both with 
the other mathematical models, using Morris Lecar and Izhikevich neurons, and with simulated Local Field 
Potentials (LFP). 



* Correspondence: alatteri@diees.unict.it 

^DIEEI - Universita di Catania, v.le A. Doria 6, Catania, Italy 

Full list of author information is available at the end of the article 

O© 201 1 Latteri et al; licensee BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons 
BIoIVIGCI Central Attribution License (http://creativecommons.Org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in 
any medium, provided the original work is properly cited. 



Latteri et al. Nonlinear Biomedical Physics 201 1, 5:2 
http://www.nonlinearbiomedphys.eom/content/5/1/2 



Page 2 of 14 



Background 

PD is a degenerative disorder of the central nervous sys- 
tem that often impairs motor skills, speech, and other 
functions [1]. It is characterized by muscle rigidity, tre- 
mor, a slowing of physical movements (bradykinesia) 
and, in extreme cases, a complete loss of physical move- 
ment (akinesia). The primary symptoms are the results 
of decreased stimulation of the motor cortex by the 
basal ganglia and other brain stem structures, tradition- 
ally considered as a consequence of the insufficient for- 
mation and action of dopamine, which is produced in 
the dopaminergic neurons of the Substantia Nigra reti- 
culata (SNr). Other symptoms may include high level 
cognitive dysfunction and subtle language problems, 
postural instability and gait disturbances. In some cases, 
it would be inaccurate to say that the cause is 
"unknown", because a small proportion is caused by 
genetic mutations. It is possible for a patient to be initi- 
ally diagnosed with PD but then to develop additional 
features, requiring revision of the diagnosis [2]. 

At present, there is no cure for PD, but medications 
or surgery can provide relief from the symptoms. The 
most widely used form of treatment is L-dopa in various 
forms. However, only 1-5% of L-DOPA enters the dopa- 
minergic neurons. The remaining L-DOPA is often 
metabolized to dopamine elsewhere, causing a wide vari- 
ety of side effects. Due to feedback inhibition, L-dopa 
results in a reduction in the endogenous formation of L- 
dopa, and so eventually becomes counterproductive [3]. 
In the 1990s the surgical ablation has been used to treat 
PD. Ablative brain surgery is the surgical lesion by 
burning or freezing or with chemical substances of brain 
tissue to treat neurological or psychological disorders. 
The thalamus was a potential target for treating tremor, 
especially in the Ventral Intermediate Nucleus (VIM) 
and Centrum Medianum-Parafascicular (CM-Pf) nuclei. 
The lesions caused by this type of surgery are irreversi- 
ble, so generally DBS surgery is considered preferable to 
lesion because it has the same effect and is adjustable 
and reversible [4-6]. Treating PD with surgery was once 
a common practice, but after the discovery of levodopa, 
surgery was restricted to only a few cases. DBS is a sur- 
gical treatment involving the implantation of a medical 
device called a brain pacemaker, which sends electrical 
impulses to specific parts of the brain. DBS in selected 
brain regions has provided remarkable therapeutic bene- 
fits for otherwise treatment-resistant movements and 
affective disorders such as chronic pain, PD or essential 
tremor and dystonia [7]. Despite the long history of 
DBS, [8], its underlying principles and mechanisms are 
still unclear. DBS directly changes brain activity in a 
controlled manner, its effects are reversible (unlike 
those of lesioning techniques). The Food and Drug 
Administration (FDA) approved DBS as a treatment for 



essential tremor in 1997, for PD in 2002 [9], and dysto- 
nia in 2003 [10]. DBS leads are placed in the brain 
according to the type of symptoms to be addressed. For 
dystonia and symptoms associated with PD (rigidity, 
bradykinesia/akinesia and tremor), the lead may be 
placed in either the Globus Pallidus or Subthalamic 
Nucleus [11]. The right side of the brain is stimulated 
to address symptoms on the left side of the body and 
vice versa. DBS does not cure PD, but it can help to 
manage some of its symptoms and subsequently to 
improve the patient's quality of life [12]. Presently, the 
procedure is used only for patients whose symptoms 
cannot be adequately controlled with medications, or 
whose medications have severe side effects [13]. Its 
direct effect on the physiology of brain cells and/or neu- 
rotransmitters is currently debated, but it is apparent 
that sending high frequency electrical impulses into spe- 
cific areas of the brain can mitigate symptoms [14] and/ 
or directly decrease the side effects induced by Parkinso- 
nian medications, [15], since it allows a sometimes huge 
decrease in medications, making the medication regime 
more tolerable. More recent neurophysiological data 
suggest that the DBS can modify also the connections 
among the cellular networks, giving rise to a holistic 
interpretation of DBS action these aspects that should 
be also considered in the future [16]. 

There are a few sites in the brain that can be targeted 
to achieve different results, so each patient must be 
assessed individually, and the particular site (or concur- 
rent sites) to be stimulated is chosen based on the parti- 
cular needs. Traditionally, the two most common sites 
are the Subthalamic Nucleus (STN) and the Globus Pal- 
lidus internus (GPi), but other sites, such as the CM-Pf 
[17] and Nucleus Tegmenti Peduncolopontini (PPTg) 
have recently shown important benefits for PD treat- 
ment [18]. A tailored DBS with the targeting of multiple 
nuclei was proposed to obtain the high clinical results 
in complex parkinsonian syndromes [18]. The main 
objective of stimulation is to re-establish desynchroniza- 
tion via a pulse train, whose parameters are selected by 
the Neurosurgeon at the main aim to decrease the dis- 
ease symptoms. The patient with PD has a high syn- 
chronization into band P (10^30 Hz) of the activities 
STN and GPi (brakinesia., akinesia., etc.), whereas the 
patient with L-Dopa into same band has high desyn- 
chronization, with an improvement of movements [19]. 
These results confirm what is expected from the Gate 
Control Theory [20], which states that the synchroniza- 
tion of neuronal activity obstructs information flow in 
brain structures, whereas, the desynchronization allows 
information flow. 

A increase of the neuronal synchronization causes a 
increase of power spectral density (PSD) of the Local 
Field Potential (LFP) taken from STN and GPi, whereas. 
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the neuronal desynchronization causes a decrease of the 
PSD. 

Such a desynchronization in the P band takes place 
also when the patient is subjected to DBS, sending a 
pulse train with frequency larger than 70 Hz [21]. The 
DBS of the STN and GPi, involves a synchronization of 
neurons in the high frequencies and a desychronization 
in the p band. These stimuli have the common goal of 
reducing the synchronized activity of specific clusters of 
neurons, in order to desynchronize the target popula- 
tion, re-establishing a normal physiological activity in a 
highly synchronized population. 

The first step toward the understanding of the 
dynamic effect of the DBS on the brain nuclei, is to 
design a model of these areas and simulate the effects 
caused by the stimulation. At this aim, after a deep 
scanning of the state of the art, in [22], a model of the 
neuronal BG was derived using a neural network made 
up of Morris-Lecar neurons, arranged in a specific 
population and showing, without any external input, 
basic oscillations mimicking the Parkinsonian tremor. 
As in [22], the positive effect of an electrical stimulation 
induced by implanted electrodes was demonstrated in 
such a population model. Therefore this model was 
initially taken into consideration as a reference. Then a 
fundamental issue, when the need of simulating large 
scale populations arises, lies in the fact that the number 
of floating point operations needed to simulate a time 
unit for a single neuron is relevant for the time needed 
to appreciate the results of the whole population. This 
led us to analyze other neuron models that could reveal 
the same characteristics (basically autonomous bursting 
oscillations) as the Morris Lecar model, but being at the 
same time much less demanding as regards the compu- 
tational burden. Therefore, after analysing, in the first 
part of this paper in the method section, the effect pro- 
duced by a population of Morris Lecar units, the same 
conditions were reproduced using a population consist- 
ing of Izhikevich neurons [23], where parameters were 
selected so as to show a resonant-like bursting charac- 
teristics, similar to the Morris Lecar model, but with a 
much less computational power demand. Of course, all 
the other characteristics, like the synaptic dynamics and 
all the other relevant parameters, were left unchanged. 

In the second part of this paper, this new reduced 
order model gave us the possibility to increase, where 
needed, the population size using roughly the same 
computation time. The positioning of the electrodes 
within the population so built revealed the expected 
effects i.e. a desynchronization within the neurons, 
which led to a spreading of the power spectral density 
over a wide frequency range. This effect even increases 
in the presence of multiple stimulating sites within the 
population. 



In this paper the population targets simulated repre- 
sented the (STN) and Globus Pallidus Externus (GPe) 
but recent researches have a new target like the PPn 
[16]. 

Methods 

The patients series and selection criteria, surgical plan- 
ning and stereotactic technique are described elsewhere 
[18]. In this paper, as also reported above, attention is 
focused on the reduced order mathematical models of 
neural structures able to show pathological PD condi- 
tions and to show the positive effect of electrical stimu- 
lation. The first mathematical model that we used 
mimics some important characteristics of the dynamical 
behaviour of a population of neurons of the STN, e.g. a 
burst and spike activity. In particular, our model dis- 
plays oscillatory activity, as observed in experimental 
investigations of brain areas, such as the basal ganglia 
and the STN, relevant for the characteristics for PD 
[24]. In this method we used the well known Morris- 
Lecar equation as a spike generator [25]. The single 
STN neuron was described by a set of four time-delayed 
differential equations, whose parameters and their role 
in the overall neuron/connection dynamics are described 
below: 

= -gcaniinfiVj - Vca) - gkWjiVj - Vfe) - gj[Vj - V/) + if""" + + I'j""'' 



dWj ^ ^ [WinfjVj) - Wj] 
dt Tuj[Vj) 

-±- = Sj[{v*-Vj{t -rj))-alfn: 
dsf- 1 

where Vj denoted the membrane potential of the jth 
neuron and Wj was an auxiliary variable Eq. 4 describes 
the evolution of the parameter ^ (the index s denotes 
synaptic), relevant for the synaptic coupling and modu- 
lating the dynamics of the current lj^[t) (see Eq. 10). 
The ^ is calculated as a function of the membrane 
potential Vy. The other parameters denote: 

C: membrane capacitance; 

gca> gk, Si' leak, Ca++, and K+ conductances through 
membranes channel; 

^ca> ^k, ♦ equilibrium potential of relevant ion 
channels; 

^2, ^3, ^4,'- tuning parameters for steady state and 
time constant. 
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The dynamics of the calcium and potassium ion chan- 
nels, responsible for the neurons spiking were controlled 
by a set of three equations: 



mnfM = 0.5 [1 + tanh{y - V1/V2}]; 
^infi^) = 0.5 [1 + tanh{y — V3/V4}]; 
Tuj[v) = l/cosh{(y - V3)/(2V4)}; 



(5) 
(6) 
(7) 



The neuron model is dimensionless, and we modu- 
lated the parameters to obtain realistic bursting patterns 
[26,27]. 

The dynamics of the neuron was controlled by exter- 
nal currents: a slowly varying current 1^'^^, whose 
dynamics is described in Eq. 3; this current had been 
proposed by Rinzel and Ermentrout [28] as a source of 
bursting, which reflects the inhibitory feedback from the 
GPe. 

The GPe neurons were excited by the STN activity Vp 
and after a time delay Xy the activity in the GPe, in turn, 
led to an inhibition of the STN neurons mediated by 
recurrent pathways [29-33]. The bursting behaviour and 
its frequency were controlled, via Eq. 3 by the Gaussian 
distributed parameters Sj and by the parameter a. The 
former induces a slightly different natural frequency 
within each neuron, whereas the latter modulates the 
inhibitory response of the GPe neurons. Background 
activity introduced by external and internal sources was 
modelled by a spatially incoherent exponentially corre- 
lated noise source with amplitude D^oise [34]. The expo- 
nentially correlated noise was calculated by using a 
second order algorithm with decay time x^o/se- These 
two last parameters contribute to the definition of 
The /"^^^^ was described by a set of two time-delayed dif- 
ferential equations: 



dt 



dXf 



— ^^lDnoise^{t)) 



dt 



, f vnoise jh 



noise y 



(8) 



(9) 



where ^(t) is Gaussian distributed, zero mean noise 
with unit standard deviation [34]. The neuronal 
dynamics displayed by the simplified model STN neu- 
rons are characterized by a bursting activity, wherein 
the simulations the bursting frequency is chosen to be 
close to frequencies usually observed in patients suffer- 
ing from parkinsonian tremor. The bursts are formed by 
a small number of spikes, six to ten per burst, which are 
controlled by the feedback from the GPe. 



The bursting pattern displayed by the network ele- 
ments is qualitatively similar to the bursting pattern 
observed in wide areas of the nervous system 
[35,36,26,27]. For illustration, we represent the neurons 
in the populations as arranged in square lattices. Two 
stimulation protocols were used: the first one using one 
electrode positioned at the lattice center, and a second 
one consisting of four stimulation electrodes, equally 
spaced within the population. 

For the numerical integration of the system we used 
the routine by Runge-Kutta method of order four. 

Simplifying the complex organization of the neuronal 
population of the STN, we mimic excitatory couplings 
between different neurons of the population. In litera- 
ture, the origin of the synchronized activity established 
between neurons in pathological conditions is not 
known, but could be induced by the modification of the 
bursting activity of the neurons. 

Also, in our model, single spikes are not able to cause 
a significant change of the spiking times of the other 
neurons, only synchronized bursts are able to induce 
such a change in the firing pattern. 

To design the connection within STN neurons we 
used a simple excitatory coupling (ij^) that does not 
take into account any input from other BG nuclei. The 
synaptic interaction is modelled as suggested by Terman 
[37,38]. 

The action potential results in an opening of the cor- 
responding ion gates, causing the income of the ij^ cur- 
rent (Eq. 10). This is composed of the local gating 
variables ^, weighted with a distance-dependent func- 
tion and multiplied with a maximal gating term and the 
potential difference corresponding to the glutamatergic 
synapses present in the STN [38-40]: 



-syn 



(0 =gsiVj - Vs) 



lly,-yfell^ 
2a} 



(10) 



Here | Ij; - j 1 1 is the distance between kth the jth and 
the neuron. The details of connections within the 
populations on which most of the current studies of 
PD focus (the BG and especially the STN) are poorly 
understood [37,40,41,38]. However, from the hippo- 
campus and the visual cortex, we know that local 
rather than global connections are developed [42,43]. 
Hence, we used a local pattern of synaptic interaction, 
cn is a normalization factor. N is the number of neu- 
rons within the population. The neurons were 
arranged on a square lattice (lattice distance dl), as 
shown in Figure 1. 

The stimulation was applied via one or four electro- 
des located within the network. The strength of the 
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Figure 1 Synchronized activity of the coupled neuron 
population. The membrane potentials of four neurons: V/ftj v 2i(t) 
V5o(t) (t), randomly selected within a regular lattice of 100 
neurons. Subscripts indicate their position in the lattice, ordered 
lexicographically. The membrane potential dynamics show the 
characteristic synchronized bursting activity, whereas in Figure 1, 
bottom side, the plot of the membrane potential on two different 
randomly selected neurons within the network is depicted, 
outlining the large amount of synchronization. The upper part of 
the diagonal shows a small side dispersion, due to the slightly 
different times of the bursting activity. Typical numerical values for 
the different variables are as follows: Vj(t)e [-0.3,0.14]; Ijsiow (0 e 
[-0.008,0.012; ^ G [0, 0.0065]. Parameters: gca= 1.0; l.Oig, = 
0.5;Vca= 1.0; v k= - 0.7; Vi= -0.5; i/? = -0.01; V2 = 0.15; V3 = 0.1; V4 = 
0.145; C= 1.0; 0 = 1.15; = -0.22; a = 0.0; tj = 10; Gaussian 
distributed sj = 2-10-3(± 2-10-5) mean (± standard deviation); as = 



dl = 0.1; Dn, 



■ 0.00001; Tnc 



: 5; Cs = 0.1; Cn= 2.0584. 



Stimulation typically decays with distance from the sti- 
mulation electrode. Which compartments of a neuron 
are activated by an electrical extracellular stimulation 
is still not clear [44,45]. The final result of the stimula- 
tion is most probably a combination of excitatory 
action directly at the soma of the neuron together with 



an activation of afferent fibres [46,31]. Therefore in 
our mathematical model we elaborated excitatory 
actions, which take place at the soma and within the 
population mediated by excitatory synapses connecting 
the STN neurons [47,40,30,41,48], and inhibitory 
actions coming from an activation of inhibitory affer- 
ent fibres [30-33]. 

Hence, the absolute value of the stimulation current 
was used to control the slowly varying current. There- 
fore, this type of stimulation mimicked the activation of 
afferent inhibitory fibres. The exact interplay between 
these different types of action remains unknown, to 
date, and is essential in determining the shape of the 
electrical pulses [49]. To cope with this challenge we 
first investigated excitatory effects. The corresponding 
stimulation term was given in Eq. (1) as XiSiI^^^^. Here 
5i = 1 is used to weight the local effectiveness of the sti- 
mulation and the step function Xi defines the onset and 
offset of the excitatory stimulation acting directly at the 
soma of the stimulated neuron. 

For the standard high-frequency DBS we used biphasic 
stimulation pulses. The first short positive pulse of 
length 0.2 ms was followed by a longer negative pulse of 
length 3.0 ms. The amplitude of the second pulse was 
adjusted such that the biphasic pulse was charged 
balanced, e.g. the total amount of applied electric charge 
was zero (Figure 2). 

Biphasic pulses are used in clinical applications to 
avoid a charge deposition in the tissue. The stimulation 
was administered through one of four electrodes posi- 
tioned in the center of the population. The spatial acti- 
vation profile of the stimulation is not known in detail 
[50] . For simplicity, we assumed a certain level of homo- 
geneity in the network and supposed that the stimula- 
tion signal uniformly and exponentially decays with 
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Figure 2 High-frequency pulse train W(t) delivered through the 
electrode. Pulses were permanently supplied at a frequency of 130 
Hz. The figure shows the typical shape of a charge balanced pulse. 
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increasing distance to the stimulation electrode. The sti- 
mulation current was given by: 



'(0 



-llly-y'll 



Wit); 



(11) 



where | \yj - y^\ \ is the distance between the k^h elec- 
trode and neuron, and is the parameter that con- 
trols the strength of the stimulation. Xi{t) (Eq. 1) 
determines the onset and offset of the excitatory stimu- 
lation. W{t) is the continuous pulse train consisting of a 
high-frequency repetition of the biphasic pulses pre- 
sented through the electrode (Figure 2). c„ is a normali- 
zation factor, which is dependent on the number of 
stimulation electrodes used and which guarantees that 
the total amount of stimulation remains independent on 
the number of electrodes used. 

To analyze the synchronization of dynamics of the 
neuronal population, we quantified the phase of the net- 
work considering the distribution of the phases. In the 
network the busting activity was the prominent 
dynamics, so we identified the bursting onset as a trig- 
ger to detect the phase displacement, calculated, for a 
single neuron j, as follows: 



Pj[t) = 27T 



tk 



tk+i - tk 



(12) 



where t g [t/^^ t/^+i], and tf^ was the onset time of the 
kfh burst of the neuron. The quantity characterizing the 
synchronization activity in oscillatory networks is given 
by: 



1 ^ 

R{t)exp[i@{t)] = — exp[i@{t)] 



(13) 



Here R is the measure of synchronization and 0 is 
mean phase. It results that 0 < R{t) < 1 for all times t; R 
= 1 corresponds to perfect in-phase synchronization, 
whereas R = 0 means complete desynchronization [5]. 

Therefore with the synchronization measure R we 
were able to reliably detect in-phase synchronization 
and desynchronization. As an example, in Figure 3 the 
value of R(t) is reported for the whole network of 100 
neurons already discussed (see Figure 1). Due to the 
random initial conditions, R starts from an initial 
(already high) value and reaches soon its maximum. 

Stimulation consists in introducing a stimulation cur- 
rent that influences Eq. 1. In this equation we remind 
that the parameter Xi controls the onset and offset of 
the excitatory stimulation acting directly at the soma of 
the stimulated neuron. Moreover W(t) is a pulse train 
with a fixed high-frequency equal to 130 Hz, consisting 
of a repetition of the biphasic pulses delivered through 
the electrode. This approach was applied in both cases 
studied, i.e. when only one electrode is implanted at the 
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Figure 3 Synchronized activity of the coupled population of 
100 neurons. (A) Bursting pattern of membrane potentials v/t); (B) 
synclironization measure R given by Eq. 13. 



centre of the network and when four electrodes are 
positioned symmetrically with respect to the lattice. 
Multi target stimulation reflects the common neurosur- 
gical protocol. However, both the techniques (single and 
multi-electrodes) have the common goal of reducing the 
synchronized activity of the target population, either by 
re-establishing a normal physiological activity in a highly 
synchronized population of neurons or by reversibly 
mimicking a tissue lesioning (DBS). 

Stimulation methods via one central electrode 

Figure 4 reports a simulation with 100 neurons orga- 
nized in regular network with one central stimulating 
electrode. 

The figure shows that, in absence of stimulation, the 
neurons bursting was strongly synchronized in-phase. 
At ^ = 4.5 s we introduced the stimulating signal: this 
caused a desynchronization for ^ = 4.5 s , leading a large 
decrease in R(t) (average value R(t) = 0.4). The pulse 
input, simulating DBS, acted so as to mask the synchro- 
nizing effects of the excitatory interconnections within 
the STN. 

Stimulation methods via four electrodes 

in the multi electrodes case, the stimulation strength (in 
terms of current in Eq. 11) decayed with distance from 
the stimulation site and each neuron received a mixture 
of the four stimulation signals. 

The DBS with four electrodes resulted in a good 
desynchronization for excitatory stimulation, indicated 
by a vanishing synchronization measure R (Figure 5(B)). 
The mean synchronization measure was R(t) = 0.15 for 
excitatory stimulation whereas, without stimulation, the 
index of synchronization was R(t) = 0.9. 
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Figure 4 Effect of a simulation of the standard high-frequency 
deep brain stimulation, via one central electrode, in terms 
synchronized activity of the coupled population of 100 
neurons. (A) Bursting pattern during tine excitatory stimulation of all 
the membrane potentials v/t). These show the characteristic 
synchronized (before stimulation) and desynchronized (after 
stimulation) bursting activity; (B) Synchronization measure R given 
by Eq. 13. In plots (A-B), the stimulation signal W(t), (see Figure 2), 
was formed by a 130 Hz pulse train composed of biphasic pulses 
(0.2 ms positive followed by 3 ms negative stimulation). The 
stimulation was supplied via one electrode situated in the central 
part of the network. Parameters: X^= 1 for t g [4.5, 9]s, excitatory 
stimulation; Xj = 0, elsewhere, 
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Figure 5 Effect of a simulation of the standard high-frequency 
deep brain stimulation, via four central electrodes, in terms 
synchronized activity of the coupled population of 100 
neurons. (A) Bursting pattern during the excitatory stimulation; the 
membrane potentials v/t) of 100 neurons are plotted, showing the 
characteristic synchronized and desynchronized bursting activity; (B) 
Synchronization measure R given by Eq. 13. In plots (A-B), the 
stimulation signal W(t), Figure 2, was formed by a 130 Hz pulse train 
composed out of biphasic pulses (0.2 ms positive followed by 3 ms 
negative stimulation). The stimulation was supplied via four 
electrodes situated in the central network. Parameters: X^ = 1 for t 
G [4.5,9]s, excitatory stimulation; X^= 0 elsewhere. DBS: Standard 
high-frequency deep brain stimulation. 

V J 



Results 

The model of the neural network that we used to simu- 
late the biological behavior, in the previous section, was 
made up of a limited number of neurons, imposed by 
the limitations of computational resources. 

Till now we modelled the neural network using Mor- 
ris Lecar neurons consisting of six differential equations: 
in these conditions, considering a reasonable computa- 
tion burden, we could analyze a number of neurons not 
larger than 100. For this reason we decided to study a 
reduced order neuron model, with less differential equa- 
tions but preserving the same characteristics. 

We used the Izhikevich model equation, in particular 
the chattering model for bursting neurons [23]. This 
model offers a reduced complexity, if compared with the 
Morris Lecar equations, retaining, at the same time, the 
most important features. In particular the neurons can 
fire with stereotypical bursts of closely spaced spikes 
[51]. 

All the other characteristics, like the synaptic 
dynamics {Ijsyn)> the noise component {Ijnoise)> the stimu- 
lation current {Ijstim) and other relevant parameters for 
the network, were left unchanged. The Izhikevich para- 
meters were suitably modified in order to allow a (time) 
comparison with the Morris Lecar model. 



The single Izhikevich membrane neuron model was 
described by the two differential equations (14 and 15). 
In dimensionless form the dynamics of the membrane 
potential Vj of the jt^ neuron, including the dynamics of 
the synaptic coupling, taken from Eq.(4), is described by 
the following set of equations: 

^ = 0.04y/ + 5vj + 140 -Uj + nif + + " + Xslf^; (14) 

^ = aihvf - (15) 
ds"- 1 

-^-"s^—^Sl-^,)-(^s^f, (16) 

v(t) = 30mVift > Tl, then \ ^ ~^ , 
^ ^ [u = u + a 

Here Uj is an auxiliary variable. The term 
0.04y^^ + 5vj + 140 - Wjwas obtained by fitting the spike 
initiation dynamics of cortical neurons so that units of 
Vj correspond to mV and units of time correspond to 
ms. We modulated the parameters so as to obtain realis- 
tic bursting patterns. The parameter a determines the 
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rate of recovery, whereas b determines the sensitivity of 
the recovery variable Uj to the membrane potential Vj . 

The parameter c determines the after-spike reset 
value, which depends on fast high threshold conduc- 
tance. Similarly, the parameter d determines the after- 
spike reset of the recovery variable, which depends on 
slow high-threshold conductance. 

The Eq. 16, was the same as Eq.4 related the system 
studied in the previous section. The Ij^oise calculated 
offline by using a second order algorithm with decay 
time Tyioise [34]. We used a numerical integration 
method based on the Eulero algorithm. 

The parameters /J^"^^ and d control the bursting beha- 
viour and its frequency. Also in this case the spikes that 
form one burst are from 6 to 10 per burst. The 
dynamics of a single neuron described by the model 
presented is depicted in Figure 6. 

To shape the connections we used the same synaptic 
current seen in the method section to connect Morris 
Lecar neurons, except that we introduced a Gaussian 
distribution into the distance among neurons to better 
emulate the realistic case. 



(17) 



Here | \yj - yj^\ \ is the distance between the /c^/^and the 
jtii neuron and is a random variable with Gaussian 
distribution, zero mean and 0.1 standard deviation; it is 
added to the distance among neurons. 
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Figure 6 Activity of one neuron described by Izhikevich 
equations. The membrane potentials y(t) of one neuron shows the 
characteristic spil<ing and bursting activity. Typical numerical values 
for the different variables are as follows. Parameters: Dnoise = 
0.00001. n = 0.001; If"'^ = 5, a = 0.02, b -- 
22, as = 0.1; ft = 0.05; Os = 0.2; cr^ = 0.02; 



0.2;c = 50, d = 0.7, T^ 
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Figure 7 Synchronized activity of the coupled population of 
neurons. The membrane potentials v/t) of 2 on 225 selected 
neurons, in one second, show the characteristic synchronized 
bursting activity. Parameters: T^ = 22, as = 0.1; ft = 0.05; Os = 0.2; Gs 
= Omig's = 02:v^ = -0.8; c/ = 2.1765; 



N is the number of neurons in the network organized 
in this random lattice (see Figure 7). 

After the definition of the reduced neuron model and 
the corresponding network, we use the same stimulation 
protocol as that applied with the Morris Lecar system 
but with a different number of neurons. 

Using numerical methods we investigated the dynami- 
cal properties of a neural population consisting of 250 
neurons. In the simulations we used the same periodic 
train pulse discussed above. In the first case, as seen 
before, we stimulated with an electrode in the middle of 
the network and then with four electrodes positioned in 
the centre of the population. A further step was to con- 
nect two identical populations to appreciate the effect of 
the stimulation in both populations, resembling the 
effect of multiple nuclei, typical of DBS. 

Also in this case the R parameter was used as the 
index of synchronization (Eq. 13). In this section we 
study the synchronization and desynchronization with 
the same protocol as in method section. So in the first 
equation with the parameter X we can determine the 
onset and offset of the excitatory stimulation, with s = 
50 and W(t) as explained above. In Figure 8 a popula- 
tion of 225 Neurons organized in random network was 
stimulated with a central electrode. 

This shows that the neuron bursting was strongly syn- 
chronized in-phase and R(t) = 1 for t g 0[3]s. The mean 
synchronization measure was R(t) = 0.6 for t g [3,6.5] s. 
In this case the network desynchronized after applying 
the stimulation signal. 

The DBS with four electrodes resulted in a good 
desynchronization for excitatory stimulation, indicated 
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Time(ms) 

Figure 8 The membrane potentials v(t) of one population of 
225 neurons organized in a random network and stimulated 
via one electrode positioned in the middle of the network. 

X = 1 for t G [3,6.5]s, excitatory stimulation; X = 0 elsewliere. 



by a vanishing synchronization measure R. In this case, 
(see Figure 9), the mean synchronization measure was R 
(t) = 0.15 for excitatory stimulation whereas without sti- 
mulation the index of synchronization was R(t) = 1. 

These results, compared with those ones reported in 
the previous section, allow us to appreciate the same 
dynamical effect of the stimulation signals also in this 
reduced order network. In summary we used: 




5000 6000 




Time(ms) 

Figure 9 The membrane potentials v(t) of 225 neurons are 
plotted: they show the characteristic synchronized and 
desynchronized bursting activity. Tine population was organized 
in random network and stimulated via four electrodes. R: 
Synchronization measure X = 1 for t g [3,6.5]s, excitatory 
stimulation; X = 0 elsewhere. 



♦ Different, reduced order model consisting of Izhi- 
kevich- Chattering; 

♦ Number of neurons equal to 225; 

♦ Random network. 

With this different system we obtained good values of 
index synchronization. For this reason we decided to 
continue the study with two coupled populations to 
study the mutual interaction effects. 

In this stimulation protocol we consider a setup, 
where a strongly synchronized neuronal population 
(population 1) acts as a pacemaker and drives another 
population (population 2), which gets synchronized 
only because of the driving. This structure resembles 
the case in which the pacemaker-like population in 
the BG and thalamus drives cortical motor 
areas which induce the peripheral shaking [22]. 
Consequently we modelled two neuronal populations 
(Figure 10): a driver (pacemaker) and a population 
(cortex) driven by the pacemaker via synaptic connec- 
tions. Within each population the coupling is local, 
respectively, whereas the coupling strengths between 
the two populations are randomized and obey a Gaus- 
sian distribution. 

To study the challenging situation of strong driving, 
we assume that the mean coupling within the driving 
population is equal to the mean coupling between the 
two populations. The excitatory coupling between the 
pacemaker and the driven system is given by ^5 = 0.4. 
Within the driven population a weak excitatory synaptic 
coupling exists, which by itself does not induce synchro- 
nization. For illustration, we represent the neurons in 
the populations as arranged in square lattices and stimu- 
lated either via four stimulation electrodes equally 
spaced within the population or with one electrode situ- 
ated in the central part of the network. 

The local separation of stimulation and recording 
sites guarantees that the feedback signal is not cor- 
rupted by stimulation artifacts. Figure 11 and 12 show 
the two populations without stimulation but coupled 
together. 




populations 



population 1 

Figure 10 Schematic plot of the stimulation setup. 
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Figure 1 1 The membrane potentials v(t) of the population 1 of 
225 neurons connected with population 2 Figure 12. The 

population was organized as a random networl< and not stimulated. 
R: Synchronization measure X = 0 at all times. 

V J 



The stimulation in population 1 causes an instanta- 
neous desynchronization of population 2, which shows a 
decrease of the synchronization measure R. 

In this case we use only one electrode to stimulate the 
coupled networks (Figure 13 and 14). 

We can see that, when stimulated, the first population 
acts as peacemaker for population 2 (driven population). 
Figure 13, resulting in a decrease of the synchronization 
index R. For the first population R(t) = 0.7 (where t g 
[3,6.5] s). For the second population. Figure 14, the value 
of R changes from an average of R(t) = 0.7 g (when t 0 
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Figure 12 The membrane potentials v(t) of the population 2 of 
225 neurons connected to population 1 (see Figurell). The 

population was organized as a random network and not stimulated. 
R: Synchronization measure X = 0 at all times. 

v J 
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Figure 13 The membrane potentials v(t) of the population 1 of 
225 neurons connected with population 2 (see Figure 14). The 

population was organized as a random network and stimulated via 
one electrode. R: Synchronization measure X = 1 for t g [3,6.5]s, 
excitatory stimulation; X = 0. 

V J 

[3]s), to R(t) = 0.5 (when t g [3,6.5] s), as a consequence 
of the coupling with population 1. 

The stimulation is realised via four electrodes located 
within the network. So we studied how the level of 
desyncronization of the driven population changes with 
a number of electrodes applied (Figure 15 and 16). 

We can see in this simulation how the increasing the 
number of electrodes causes an increasing in the desyn- 
chronization effect. Stimulating the first population (Fig- 
ure 15) acts as peacemaker for the driven population 2 
(Figure 16). 




0 1000 2000 3000 4000 5000 6000 



Time(ms) 

Figure 14 The membrane potentials v(t) of the population 2 of 
225 neurons connected with population 1 (see Figure 13). The 

population was organized as a random network and not stimulated. 

R: Synchronization measure. 
J 
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Time(ms) 

Figure 15 The membrane potentials v(t) of the population 1 of 
225 neurons connected with population 2, Figure 16. The 

population was organized as a random network and stimulated via 

four electrodes. R: Synchronization measure X = 1 for t [3,6.5]s, 

excitatory stimulation; X = 0 elsewhere. 
\ / 



For the first population we have R(t) = 0.2 where t g 
[3,6.5] s, moreover for the second population the value of 
R change from R(t) = 0.5 with t g 0 [3]s, to R(t) = 0.3 
with t G [3,6.5] s, given by coupling with population 1. 

Discussion 

The macroscopic effect of stimulation is represented by 
a modification of the Spectrograms with respect to the 
non-stimulated case. To appreciate the introduced 
model and its response to a stimulating pulse train, we 
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Time(ms) 

Figure 16 The membrane potentials v(t) of the population 2 of 
225 neurons connected with first population ^, Figure 15. The 

population was organized as a random network and no stimulated. 
R: Synchronization measure. 

v J 



calculate the LPF for our network model and translate 
the single burst activity into a macroscopic field mea- 
sure. The LFP was calculated by an electrode positioned 
in the center of the neuronal population and was given 
by: 

no = ^E^; (18) 

where rj is the distance between neuron ; and the 
recording electrode [50]. All the current components, 
except for the stimulation current, contributed to the 
ionic current Ij{t)(Eq. 1). Therefore Ij{t) is composed of 
jconst^ jconst^ jslow^ fw ^^^^^ .^^.^ currents respon- 

sible for the basic spiking features of the jth STN model 
neurons (i.e., the right hand side of Eq. 1). Rg = 1 is the 
extracellular resistivity that was assumed to be 
homogeneous. 

The dynamics of the non stimulated network for the 
Moris Lecar model, in terms of LFP is reported (through 
the spectrogram in dB) in Figure 17. From this plot we 
can appreciate the presence of a high synchronization at 
low frequencies: this is the key characteristic of a patho- 
logical network, according to the neuronal gate theory. 
It is to be outlined that we conducted the study consid- 
ering a very low number of neurons with respect to the 
actual case. Moreover relevant parameters, like the noise 
level and sources and the details of the network topol- 
ogy, in the in vivo case, were mostly unknown. Under 
this perspective, the results obtained can be considered 
relevant, above all in view of the effect of stimulation. 

In the following the effect of stimulation is considered. 
We apply the traditional stimulation via four electrodes 
in the same 100 neuron population. 




Time s 



Figure 17 Spectrogram in dB. LFP of a not stimulated population 

of 100 Morris Lecar neurons. 
) 
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When we apply the pulse train at 130 Hz the network 
desynchronizes particularly in the low frequency range, 
leading to a certain spreading of the spectral intensity, 
as we can see in Figure 18 for the Morris Lecar model 
The results shown in Figure 18 can be compared with 
the literature [21] when the patient is treated with L- 
Dopa. These results were mostly important for our 
study: both in the un-stimulated and in the stimulated 
case, we can appreciate the same effect of stimulation, 
both in simulation and experimentally. 

Referring to one population, in order to compare the 
dynamics of the un-stimulated and the stimulated net- 
work of the Izhijevich and Morris Lecar network, we 
studied how the application of a stimulation signal 
caused a change in the LFP of a population. As seen 
before, we will appreciate a desynchonization of the LFP 
through the spectrogram in dB of Izhikevich network. 
From Eq.l8 we calculated the LFP of the Izhikevich 
population. 

Now we can compare this LFP simulated signal of 
Izhikevich neuron population with the LFP simulated 
signal of the Morris Lecar population (Figure 17, Figure 
19), both not stimulated. 

From these plots we can see that there is a high syn- 
chronization at low frequency and this is a characteris- 
tic, mimicking the pathological network. This behavior 
is clearly visible also in the Izhikevich neuron popula- 
tion: between the [5-20] Hz range there is a high 
activity. 

Subsequently we applied the traditional stimulation via 
four electrodes in a population (see Figure 20). 

When we apply the pulse train at 130 Hz, the network 
desynchronizes at the low frequencies and we notice a 
spreading of the power density. The same effect is seen 
in the literature, when the patient is treated with L- 
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Figure 18 Spectrogram in dB. LFP of a stimulated population for 
all times with train pulse at 130 Hz of 100 Morris Lecar neurons. 

v J 
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Figure 19 Spectrogram in db. LFP of a not stimulated population 

of 225 Izhikevich neurons. 

^ J 



Dopa [21], and in Figure 18, referring to the Morris 
Lecar network. If we compare the Figure 17 with Figure 
19 and 18 with Figure 20, we can see that we obtained a 
same effects both in the stimulated and in the not sti- 
mulated case. 

Conclusions 

In this work some models for studying the effect of sti- 
mulation in neural populations mimicking the dynamics 
met in the relevant brain neural tissues of PD patients 
are analysed. A new structure, built using a reduced 
order neuron model (the Izhikevich one), computation- 
ally much lighter, was demonstrated to lead to the same 
results as other neural networks already introduced in 
the literature. This gives the opportunity to more 
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Figure 20 Spectrogram in db. LFP of a stimulated population for 
all times with train pulse at 130 Hz of 225 Izhikevich neurons High- 
frequency pulse train W(t) delivered through the electrode. 
J 
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efficiently study the effect of stimulation in large scale 
neural networks. Results were compared with the other 
mathematical models, i.e. Morris Lecar neurons, analyz- 
ing both parameters related to the neural level (the syn- 
chronization effect) and indexes related to the 
macroscopic level (LFP). The results agree in terms of 
the positive effect of the stimulation in terms of power 
spectral density. We consider these results very promis- 
ing: they open the way to the possibility of simulating 
large scale networks in pathological conditions at the 
aim to design new control strategies for the PD effect 
mitigation. 
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